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Abstract 

In this paper we present an improved lattice Boltzmann model for compressible 
Navier-Stokes system with high Mach number. The model is composed of three 
components: (i) the discrete-velocity-model by Watari and Tsutahara [Phys Rev E 
67,036306(2003)], (ii) a modified Lax-Wendroff finite difference scheme where rea- 
sonable dissipation and dispersion are naturally included, (iii) artificial viscosity. 
The improved model is convenient to compromise the high accuracy and stabil- 
ity. The included dispersion term can effectively reduce the numerical oscillation at 
discontinuity. The added artificial viscosity helps the scheme to satisfy the von Neu- 
mann stability condition. Shock tubes and shock reflections are used to validate the 
new scheme. In our numerical tests the Mach numbers are successfully increased up 
to 20 or higher. The flexibility of the new model makes it suitable for tracking shock 
waves with high accuracy and for investigating nonlinear nonequilibrium complex 
systems. 

Key words: Lattice Boltzmann Method, Compressible Flows, von Neumann 
Analysis 
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1 Introduction 



Recently, the lattice Boltzmann(LB) method got substantial progress and has 
been regarded as a promising alternative for simulating many complex phe- 
nomena in various fields[l]. Unlike the macroscopic computational fluid dy- 
namics or the microscopic molecular dynamics, the LB uses a mesoscopic 



Preprint submitted to Elsevier Science 



4 February 2008 



discrete Boltzmann equation to describe the fluid system. Because of its in- 
trinsic kinetic nature, the LB contains more physical connotation than Navier- 
Stokes or Euler equations based on the continuum hypothesis[2]. From the 
Chapmann-Enskog analysis, the latter can be derived from the former under 
the hydrodynamic limit. 

Although having achieved great success in simulating incompressible fluids, 
the application of LB to high-speed compressible flows still needs substantial 
effort. High-speed compressible flows are ubiquitous in various flelds, such as 
explosion physics, aeronautics and so on [3]. Simulation of the compressible 
Navier-Stokes system, especially for the those containing shock waves or con- 
tact discontinuities, is an interesting and challenging work. Along the line, ex- 
tensive efforts have been made in the past years. Alexander, et al[4] presented 
a model where the sound speed is selectable; Yan, et al[5] proposed a com- 
pressible LB model with three-speed-three-energy-level for the Euler system; 
Yu and Zhao [6] composed a model for compressible flows by introducing an at- 
tractive force to soften sound speed; Sun [7, 8, 9] contributed a locally adaptive 
semi-discrete LB model, where the set of particle speed is chosen according 
to the local fluid velocity and internal energy so that the fluid velocity is no 
longer hmited by the particle speed set. In the development of LB for Navier- 
Stokes systems, another way is referred to as the finite difference lattice Boltz- 
mann method (FDLBM) [10,11,12]. The one by Watari-Tsutahara (WT) is typ- 
ical [11]. The same idea was then extended to binary compressible fiows[12]. 
FDLBM[11,12] breaks the binding of discretizations of space and time and 
makes the particle speeds more flexible. But similar to previous LB models, 
the numerical stability problem remains one of the few blocks for its practical 
simulation to high Mach number flows. The stability problem of LB has been 
addressed and attempted for some years [13,14,15,16,17,18,19,20,21,22,23] . 
Among them, the entropic LB method[16,17] tries to make the scheme to 
follow the if-theorem; The FIX-UP method[16,18] is based on the standard 
BGK scheme, uses a third order equilibrium distribution function and a self- 
adapting updating parameter to avoid negativeness of the mass distribution 
function. Flux limiter techniques are used to enhance the stability of FDLB 
by Sofonea, et al[19]. Adding minimal dissipation locally to improve stabil- 
ity is also suggested by Brownlee, et al[20], but there such an approach is 
not explicitly discussed. All the above mentioned attempts are for low Mach 
number flows. In this paper we present a new LB scheme for high-speed com- 
pressible flows which is composed of three components, (i) the original DVM 
by WT, (ii) an Modified Lax-Wendroff (MLW) finite difference scheme where 
reasonable dissipation and dispersion are naturally included, (iii) additional 
artificial viscosity. With the new scheme, high speed compressible flows with 
strong shocks can be successfully simulated. 

This paper is organized as follows. In section 2 the original DVM by WT is 
briefly reviewed. An alternative FD scheme combined with artiflcial viscosity 
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Fig. 1. Sketch of the discrete-velocity-model used in the present paper. 

is introduced in section 3. The von Neumann stability analysis is performed 
in section 4, from which solutions to improve the numerical stability can be 
found. Several benchmark tests are used to validate the proposed scheme in 
section 5. Section 6 presents the concluding remarks. 



2 Outline of the DVM by Watari-Tsutahara 



DVM of WT can be write as: 

%Tl ITT 

vo = 0, Vki = Vk[cos{—),sm{—)],i = 1, 2. ..8, (1) 

where subscript k indicates the A;-th group of the particle velocities whose 
speed is and i indicates the direction of particle's speed. A sketch of the 
DVM is referred to Fig.l. It's easy to prove that this DVM at least up to 
seventh rank isotropy. The evolution of the distribution function f^i with the 
Bhatanger-Gross-Krook approximation [24] reads, 

^ + ->^^■^ = -l[^^-a, (2) 

where is the discrete version of the local equilibrium distribution function; 
r is the spatial coordinate; r is the relaxation time; the local particle density 
n, hydrodynamic velocity u and temperature T are defined by 
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where P and Cint are the local pressure and internal energy. This model is 
designed to recover the following Navier-Stokes equations 
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in the hydrodynamic limit, where p ,k are viscosity coefficient and heat con- 
ductivity coefficient, having the following relations with pressure P and relax- 
ation time T : 

^ = Pr, K = 2p. (9) 
The equilibrium distribution function f^^ is calculated in the following way. 
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3 Modified Lax-Wendroff scheme and cirtificial viscosity 



To simplify the discussion, we work on the general Cartesian coordinate. The 

combination of the above DVM and the general FD scheme with first-order 
forward in time and second-order upwinding in space composes the original 
FDLB model by WT. It has been validated via the Couette flow, small Mach 
number Riemann problems. When the Mach number M exceeds 1, the original 
LB model is not stable. The DVM is derived independent of Mach number. 
Therefore, we resort to the discretization of the left-hand side of Eq. (2) to 
make an accurate and stable LB scheme. Here we investigate a mixed scheme 
which is composed of a modified Lax-Wendroff[25] and an artificial viscocity. 



As we know, the original Lax-Wendroff (LW) scheme is very dissipative and 
has a strong "smoothing effect". Obviously, it is not favorable when needing 
capture shocks in the system. To compromise the accuracy and stability, we 
add a dispersion term and the artificial viscosity to the right-hand side of Eq. 
(2) before discretization so that we have 
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where 



(14) 

(15) 
(16) 

Pal+l + 2Pq;/ + PaI-1 

plays a role of the switching function, A is a coefficient controlling the ampli- 
tude of the artificial viscosity. Using the Lax-Wendroff to the left-hand side 
and central difference to the right-hand side of Eq. (14) results in the following 
LB equation. 
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where the third suffixes / — 1, /, / -|- 1 indicate the mesh nodes in a; or y 
direction. The positions of terms 3 and 4 in the right-hand side of Eq. (17) 



5 



have been exchanged. It is clear that the first fine corresponds to the general 
LB equation with the central difference in space; compared with the central 
difference, the Lax-Wendroff contributes an extra line II; lines III and IV show 
the added dispersion term and artificial viscosity. 



4 von Neumann Stability Analysis 



We analysis the numerical stability of the FDLBM by means of von Neumann 
stabihty analysis [2 1,22]. In the analysis solution of finite-difference equation 
is written as the familiar Fourier series, and the numerical stability is evalu- 
ated by the magnitude of eigenvalues of an amplification matrix. The small 
perturbation Afki is defined as: fki{r,t) = Afki{T,t) + /^j, where f^^ is the 
global equilibrium distribution function which is a constant, depends only on 
the mean density, velocity and temperature. Prom equation(14) we can obtain 

At ""''"^^ " ~r f ~ ^"'^ 

vU^ - cUArl d^Ah, " Ari d' Ah, 

+ 6 + K\ (1 - \^-\)2Ai^^- 

The perturbation part Afki{ra,t) may be written as series of complex expo- 
nents, Afki{ra-,t) = Flj^exp^ikara), where F^j is an amplitude at grid point 
and time t, i is an imaginary unit, and ka is the wave number of sine wave 
in the domain with the highest resolution l/Ava- Substituting this expansion 
into the equation (18), we obtain F^^^^ = GijFl- , where Gij is a matrix being 
used to assess amplification rate of F^^ per time step At. If the maximum of 
the eigenvalues of the amplification matrix satisfies the condition, max|a;| < 1, 
for all wave numbers, the FD scheme is surely stable, where uj is the eigenvalue 
of the amplification matrix. This is the von Neumann condition for stability. 
The amplification matrix Gij can be written as following. 



-'12 

dfl^dfld^^dfldl^^dfldu^^ 

dfkj dp dfkj dT dfkj du^ dfkj ' 
Several researchers have analyzed the stability of the incompressible LB models[21,26,27], 
it is found that there is not a single wave-number being always the most unsta- 
ble. For the 2D DVM by WT, Gij is a matrix with 33 x 33 elements. Moreover, 
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Fig. 2. (Color online) Stability analysis for four conditions. The macroscopic vari- 
ables are set as {p,ui,U2,T) = (2.0,30.0,0.0,1.0), the other constants are set as 
{vo, ?;2, ?^3, ?^4) = (0.0, 1.0, 1.92, 2.99, 4.49), A = 0.5, Ax = Ay = 3 x IQ-^At = r 
= 10-5. 

every element is related to the macroscopical variables (density, temperature, 
velocities), discrete velocities and other constants, so it is difficult to analyze 
with explicit expressions. We resort to using the software, Mathematica-5 to 
conduct a series of quantitative analysis. Now we show some numerical results 
of von Neumann analysis by Mathematica-5. The results will be shown by 
figures with curves for the maximum eigenvalue |ci;| max of Gij versus kAx. 

Figure 2 shows a comparison of the fours different cases, (i) LB with only 
the central-difference to the convection term, i.e. LB scheme based only on 
the ffist line of Eq. (17) (see the black line with squares); (ii) LB with Lax- 
Wendroff, i.e., LB scheme based on the ffist two lines of Eq. (17) (see the red 
line with circles); (iii) LB scheme based on lines 1 and 3 of Eq. (17) (see the 
green line with triangles); (iv) LB scheme based on the whole of Eq. (17) (see 
the blue line with triangles down). Here the macroscopic variables are cho- 
sen as {p,Ui,U2,T) = (2.0,30.0,0.0,1.0), and the remaining parameters are 
(t;o,^^i,^^2,^^3,?^4) = (0.0,1.0,1.92,2.99,4.49), A = 0.5, Ax = Ay = 3 x 10'^, 
At = T = 10~^. It is clear that the artificial viscosity term can significantly 
decrease the maximum eigenvalue leu] max from being larger than to be smaller 
than 1 for appropriately given time step. Moreover, it is worthy to note that 
dissipation term in line 2 of Eq. (17) favors and dispersion term in line 3 
disfavors the stability to some extent. Numerical experiments show that the 
dispersion term may effectively reduce numerical oscillations near discontinu- 
ity and improves the accuracy (see Fig. 3 for an example). 

Figure 4 shows the effects of various artificial viscosities to the stability. Fig. (a) 
shows the cases with A = 1.0, 0.5, 0.1, and 0.05. Fig.(b) shows the cases with 
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Fig. 3. (Color online) Effects of dispersion to simulation. The macroscopic vari- 
ables at two sides of a shock tube are set as {p,ui,U2,T)\l = (10.0,10.0,0.0,5.0), 
{p,ui,U2,T)\ji = (10.0,0.0,0.0,5.0), A = 2.0. The other constants and macroscopic 
variables are the same as in Fig. 2. 




Fig. 4. (Color online) Effects of various artificial viscosities to the numerical stability. 
Fig. (a) shows the cases with A = 1.0, 0.5, 0.1, 0.05. Fig. (b) shows the cases with 
A = 1.0,1.8, 2.0, 3.0. The other constants and macroscopic variables are unchanged. 

A = 1.0, 1.8,2.0, and 3.0. The other constants and macroscopic variables are 
unchanged. From this figure we can see some relevance: Strength of artificial 
viscosity has a large impact on the stability. LB works only within a certain 
range of artificial viscosity. In practical simulations, we generally take the 
smaller viscosity in favor of the accuracy. 

Since the density p can be normalized to 1, we then need only investigate the 
effects of the other two physical quantities, temperature T and flow velocity 
u. Figure 5 shows four cases with T = 1,T = 5,T = 15 and T = 25. Here 
Ml = 5, M2 = and A = 0. When other parameters are fixed, the numerical 
stability becomes better with the increasing of temperature. This can also be 
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Fig. 5. (Color online) Influence of temperature T to numerical stability, ui = 5, 
U2 = and A = 0. The other physical quantities and model parameters are un- 
changed. 
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Fig. 6. (Color online) Influence of velocity u to numerical stability. The value of 
ui is altered from zero to 30 and U2 = 0. Here A = 0, the other constants and 
macroscopic variables are unchanged. 

understood that higher temperature corresponds to higher sound speed and 
lower Mach number. 

Figure 6 shows cases with difference flow velocities. The value of Ui is altered 
from zero to 30 and U2 = 0. Here A = 0, the other constants and macroscopic 
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variables are unchanged. This figure clearly shows that the higher the Mach 
number, the larger the maximum eigenvalue, which answers why the numerical 
stability becomes worse with the increasing of Mach number of the fluid. 



5 Numerical tests and analysis 



In this section two kinds of typical benchmarks are used to validate the newly 
proposed scheme. The flrst kind is the Riemann problem[28]. The second one 
is the problem of shock reflection [29]. 



5.1 Riemann problems [28] 



Here the two-dimensional model is used to solve the one-dimensional Riemann 
problem. The initial macroscopic variables at the two sides are denoted by (p, 
Ml, M2, T)|l and (p, mi, ^2, T)\r, respectively. 



5.1.1 Sod's shock tube 

For the problem considered, the initial condition is described by 

'(p,iii,ii2,r)|L = (i.o,o.o,o.o,i.o) 

(21) 

(p,mi,M2,T)|r= (0.125,0.0,0.0,0.8) 

Figure 7 shows the computed density, pressure, velocity, temperature proflles 
at t = 0.2, where the circles are simulation results and solid lines with squares 
are analytical solutions. The size of grid is Ax = Ay = 10^^, time step 
At — 10~^, and r = 10~^, A = 2. The two sets of results have a satisfying 
agreement. 



5.1.2 Lax's shock tube 

The initial condition of this problem reads 

(p,Mi,M2,T)U = (0.445,0.698,0.0,7.928) ^^^^ 
(p, iii, 1^2, r)|fl= (0.5,0.0,0.0,1.142) 

Figure 8 shows the results at t = 0.2, where the circles are simulation results 
and solid lines with squares correspond to exact solutions. The parameters are 
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Fig. 7. (Color online) Comparison of numerical and theoretical results for the Sod 
shock tube, where t=0.2. Solid lines with squares are for exact solutions and circles 
are for simulation results. 
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Fig. 8. (Color online) Comparison of numerical and theoretical results for the Lax 
shock tube at t = 0.2. Solid lines with squares are for exact solution and circles are 
for simulation results. 



set to be Ax = Ay = 10-^ At = IQ-^, t = 10"^, and A = 1. We also find a 
good agreement between the two sets of results. 
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Fig. 9. (Color online) Comparison of numerical and theoretical results for the Sjo- 
green problem at t = 0.018. Solid lines with squares are for exact solutions, and 
circles are for simulation results. 

5.1.3 Sjogreen's problem 

The initial condition of this problem is 



(p,Mi,M2,r)U = (1.0,-2.0,0.0,0.4) 
(p,Mi,U2,T)|jj = (1.0,2.0,0.0,0.4) 



(23) 



The numerical and exact solutions at t = 0.018 are shown in Fig. 9, where 
Ax = Ay = 3 X 10-^ At = t = 10"^ and A = 1.8. The exact solution of 
this problem consists of two strong rarefaction waves and a weak constant 
contact discontinuity. Pressure near the contact discontinuity is very small, 
which brings certain difficulties to simulation. Temperature and density cal- 
culated by many schemes are negative. However, the improved model ensures 
the positivity of them. Successful simulation of this problem proves that the 
improved model is apphcable to the low-density, low-temperature flow simu- 
lations. 



5.1.4 Colella's explosion wave problem 

The initial condition of this test can be write as 

{p,UuU2,T)\l = (1.0,0.0,0.0,1000.0) 
(p,mi,M2,T)|h = (1.0,0.0,0.0,0.01) 



(24) 
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Fig. 10. (Color online) Comparison of numerical and theoretical results for the 
Colella explosion wave problem at t=0.05. Solid lines with squares are for exact 
solutions and circles are for simulation results. 



This is generally regarded as a difficult test. The exact solution contains a 
leftwards rarefaction wave, a contact discontinuity and a strong shock. It is 
generally used to check the robustness and accuracy. Figure 10 gives compar- 
ison of the numerical and theoretical results at t = 0.05. Here A = 20, other 
parameters are same as in the Sjogreen test. Successful simulation of this 
test proves the improved model is applicable to flows with very high ratios of 
temperature and pressure. 



5.1.5 Collision of two strong shocks 
This test with the following initial data: 

{p,Ui,U2,T)\l = (5.99924,19.5975,0.0,76.8254) 
{p,Ui,U2,T)\r = (5.99242,-6.19633,0.0,7.69222) 

This is also a difficult test. Exact solution contains a leftwards shock, a right 
contact discontinuity and shock which spreading to right side. And the left- 
shock spreads to right very slowly, which brings additional difficulties to the 
numerical method. Fig. 11 gives a comparison of the numerical and theoretical 
results at t = 0.12. Parameters used in this test are Ax = Ay = 2 x 10~^, 
At = T = 10~^ and A = 1. The good agreement between the two sets of results 
shows again the robustness of the improved model. 
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Fig. 11. (Color online) Comparison of numerical and theoretical results for collision 
of two strong shocks at t = 0.12. Solid lines with symbols are for exact solutions 
and symbols are for simulation results. 

5.2 Shock reflections [29] 



We will present two gas dynamics simulations. Both are done on rectangular 
grid. The first is to recover a steady regular shock reflection. The second is 
the double Mach reflection of a shock off an oblique surface. This example is 
used in Ref . [30] as a benchmark test for comparing the performance of various 
difference methods on problem involving strong shocks. 



5.2.1 Steady regular shock reflection 

In the first test problem, the incoming shock wave with Mach number 20 has 
an angle of 30° to the wall. The computational domain is a rectangle with 
length 0.9 and height 0.3. This domain is divided into a 300 x 100 rectangular 
grid with Ax = Ay = 0.003. The boundary conditions are composed of a re- 
flecting surface along the bottom boundary, supersonic outflow along the right 
boundary , and Dirichlet conditions on the left and top boundary conditions, 
given by 

(p, Ml, U2, T)|o, t = (1.0, 20.0, 0.0, 0.5), 

(26) 

^ (p, Ml, U2, T) I,, 0.3, t = (ff , 16.7, - 5.71578, 22.61). 

Figure 12 shows contours of density at t = 0.2. The clear shock reflection on 
the wall agrees well with the exact solution. 
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50 100 150 200 250 300 

Fig. 12. (Color online) Density contour of steady regular shock reflection on a wall 
at t=0.2. From black to white, the density increases. 

5.2.2 Double Mach reflection 

In this test, we considered an unsteady shock reflection. The initial pressure 
ratio here is high. A planar shock is incident towards an oblique surface with a 
30° angle to the direction of propagation of the shock. A uniform mesh size of 
500 X 200 is used for the numerical simulation. The conditions for both sides 
are: 

f (^,13.3 cos 30°, -13.3 sin 30°, 89.2775), if y > h{x,0) 

{p,Ui,U2,T) \x,yfi= { 

[(2.0,0.0,0.0,0.5) ify<h{x,0) 

(27) 

where h{x, t) = \/3{x — 80Ax) — 4:0t. The reflecting wall lines along the bottom 
of the problem domain, beginning at x = 0.08. The shock makes a 60° angle 
with the X axis and extends to the top of the problem domain ai y = 0.2. At 
the top boundary, the physical quantities are assigned the same values as on 
the left side for x < g{t) and are assigned the same values on the right side, 
where g{t) = 80 Ax + v^/3(0.2 + 40t). The computed density, temperature 
and flow velocity along the x-direction are shown in Fig. 13, where complex 
characteristics, such as obhque shocks and triple points, are well captured. 



6 Conclusions and discussions 

A lattice Boltzmann model to the high-speed compressible Navier-Stokes sys- 
tem is presented. The new LB model is composed of the following components: 
the original DVM by Watari and Tsutahara, a modifled Lax-Wendroff scheme 
and an additional artiflcial viscosity. Compared with the central difference 
scheme, the Lax-Wendroff contributes a dissipation term which is in favor of 
the numerical stability, even though it is generally still not enough for high- 
speed flows. The introducing of the third-order dispersion term helps to elim- 
inate some unphysical oscillations at discontinuity. The additional artiflcial 
viscosity compensates the insufficiency of the above-mentioned dissipation so 
that the LB simulation can continue smoothly. The adding of the dispersion 
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Fig. 13. (Color online) Contours of density (top), temperature (center), and ui 
(bottom) of the double Mach reflection problem at the 750th iteration step. The 
units of the x- and y- axes are both 0.001. 

and artificial viscosity terms should survive the dilemma of stability versus ac- 
curacy. In other words, they should be minimal but make the evolution satisfy 
the von Neumann stability condition. Due to the complexity, the analysis re- 
sorts to the software, Mathematica-5, and only some typical results are shown 
by figures. 

Typical benchmark tests are used to validate the proposed scheme. Riemann 
problems, including the Sod, Lax, Sjogreen, Colella explosion wave, collision 
of two strong shocks, show good accuracy and numerical stability of the new 
scheme, even though they are generally difficult to resolve by traditional com- 
putational fluid dynamics[28,29,30,31]. Regular and double Mach shock reflec- 
tions are successfully recovered. These simulations show that the improved LB 
model may be used to investigate some long-standing problems, such as the 
transitions between regular and Mach reflections. By incorporating an appro- 
priate equation of state, or equivalently, a free energy functional, or an external 
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force, the present model may be used to simulate liquid-vapor transition and 
relevant flow behavior. Future work includes a more complete description of 
the problem on numerical accuracy versus stability and thermal lattice Boltz- 
mann model for multi-phase flows. 
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